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Abstract. For a dielectric solid surrounded by an electrolyte and positioned 
inside an externally biased parallel-plate capacitor, we study numerically how the 
resulting induced-charge electro-osmotic (ICEO) flow depends on the topology 
and shape of the dielectric solid. In particular, we extend existing conventional 
electrokinetic models with an artificial design field to describe the transition 
from the liquid electrolyte to the solid dielectric. Using this design field, we 
have succeeded in applying the method of topology optimization to find system 
geometries with non-trivial topologies that maximize the net induced electro- 
osmotic flow rate through the electrolytic capacitor in the direction parallel to 
the capacitor plates. Once found, the performance of the topology optimized 
geometries has been validated by transferring them to conventional electrokinetic 
models not relying on the artificial design field. Our results show the importance 
of the topology and shape of the dielectric solid in ICEO systems and point to 
new designs of ICEO micropumps with significantly improved performance. 
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1. Introduction 

Induced-charge electro-osmotic (ICEO) flow is generated when an external electric field 
polarizes a solid object placed in an electrolytic solution [HE]- Initially, the object 
acquires a position-dependent potential difference £ relative to the bulk electrolyte. 
However, this potential is screened out by the counter ions in the electrolyte by the 
formation of an electric double layer of width Ad at the surface of the object. The 
ions in the diffusive part of the double layer are then electromigrating in the resulting 
external electric field, and by viscous forces they drag the fluid along. At the outer 
surface of the double layer a resulting effective slip velocity is thus established. For a 
review of ICEO see Squires and Bazant [3] . 

The ICEO effect may be utilized in microfluidic devices for fluid manipulation, 
as proposed in 2004 by Bazant and squires [1] . Theoretically, various simple dielectric 
shapes have been analyzed analytically for their ability to pump and mix liquids [31 14] . 
Experimentally ICEO was observed and the basic model validated against particle 
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image velocimetry in 2005 2|, and later it has been used in a microfluidic mixer, 
where a number of triangular shapes act as passive mixers [S]. However, no studies 
have been carried out concerning the impact of topology changes of the dielectric 
shapes on the mixing or pumping efficiency. In this work we focus on the application 
of topology optimization to ICEO systems. With this method it is possible to optimize 
the dielectric shapes for many purposes, such as mixing and pumping efficiency. 

Our model system consists of two externally biased, parallel capacitor plates 
confining an electrolyte. A dielectric solid is shaped and positioned in the electrolyte, 
and the external bias induces ICEO flow at the dielectric surfaces. In this work 
we focus on optimizing the topology and shape of the dielectric solid to generate 
the maximal flow perpendicular to the external applied electric field. This example 
of establishing an optimized ICEO micropump serves as demonstration of the 
implemented topology optimization method. 

Following the method of Borrvall and Petersson [Hj and the implementation by 
Olesen, Okkels and Bruus [7] of topology optimization in microfluidic systems we 
introduce an artificial design field 7(1-) in the governing equations. The design field 
varies continuously from zero to unity, and it defines to which degree a point in 
the design domain is occupied by dielectric solid or electrolytic fluid. Here, 7 = 
is the limit of pure solid and 7 = 1 is the limit of pure fluid, while intermediate 
values of 7 represent a mixture of solid and fluid. In this way, the discrete problem 
of placing and shaping the dielectric solid in the electrolytic fluid is converted into 
a continuous problem, where the sharp borders between solid and electrolyte are 
replaced by continuous transitions throughout the design domain. In some sense 
one can think of the solid/fluid mixture as a sort of ion-exchange membrane in the 
form of a sponge with varying permeability. This continuum formulation allows for 
an efficient gradient-based optimization of the problem. 

In one important aspect our system differs from other systems previously studied 
by topology optimization: induced-charge electro-osmosis is a boundary effect relying 
on the polarization and screening charges in a nanometer-sized region around the 
solid/fluid interface. Previously studied systems have all been relying on bulk 
properties such as the distribution of solids in mechanical stress analysis [5] , photonic 
band gap structures in optical wave guides [9] , and acoustic noise reduction [10] , or on 
the distribution of solids and liquids in viscous channel networks [HEIHT] and chemical 
microreactors j!2| . In our case, as for most other applications of topology optimization, 
no mathematical proof exists that the topology optimization routine indeed will result 
in an optimized structure. Moreover, since the boundary effects of our problem result 
in a numerical analysis which is very sensitive on the initial conditions, on the meshing, 
and on the specific form of the design field, we take the more pragmatic approach of 
finding non-trivial geometries by use of topology optimization, and then validate the 
optimality by transferring the geometries to conventional electrokinetic models not 
relying on the artificial design field. 

2. Model system 

We consider a parallel-plate capacitor externally biased with a harmonic oscillating 
voltage difference A<p = 2Vb cos(wi) and enclosing an electrolyte and a dielectric solid. 
The two capacitor plates are positioned at z = ±if and periodic boundary conditions 
are applied at the open ends at x = ±L/2. The resulting bound domain of size 
L x 2H in the rrz-plane is shown in Fig. [T] The system is assumed to be unbounded 
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Figure 1. (a) A sketch of the rectangular L X 2H cross-section of the electrolytic 
capacitor in the xz-plane. The external voltage tj>i and <j)2 is applied to the 
two infinite parallel-plate electrodes (thick black lines) at z = ±if . The voltage 
difference 01 — 02 induces an ICEO flow around the un-biased dielectric solid (dark 
gray) shaped by the topology optimization routine limited to the rectangular lx2h 
design domain (light gray). The dielectric solid is surrounded by pure electrolyte 
(light gray and white). Periodic boundary conditions are applied ar the vertical 
edges (dotted lines), (b) The dimensionless electric permittivity e as a function 
of the design variable 7. (c) Zoom-in on the rapid convergence of £(7) towards 
Efluid = 1 for 7 approaching unity after passing the value 7 cu t — A7 ~ 0.98. 



and translational invariant in the perpendicular y-direction. The topology and shape 
of the dielectric solid is determined by the numerical optimization routine acting only 
within a smaller rectangular, central design domain of size I x 2h. The remaining 
domain outside this area is filled with pure electrolyte. Double layers, or Debye 
screening layers, are formed in the electrolyte around each piece of dielectric solid 
to screen out the induced polarization charges. The pull from the external electric 
field on these screening layers in the design domain drives an ICEO flow in the entire 
domain. 

If the dielectric solid is symmetric around the x-axis, the anti-symmetry of the 
applied external bias voltage ensures that the resulting electric potential is anti- 
symmetric and the velocity and pressure fields symmetric around the center plane 
z = 0. This symmetry appears in most of the cases studied in this paper, and when 
present it is exploited to obtain a significant decrease in memory requirements of the 
numerical calculations. 

The specific goal of our analysis is to determine the topology and shape of the 
dielectric solid such that a maximal flow rate Q is obtained parallel to the i-axis, i.e. 
perpendicular to the direction of external potential field gradient. 

3. Governing equations 

We follow the conventional continuum approach to the electrokinetic modeling of the 
electrolytic capacitor [3J- For simplicity we consider a symmetric, binary electrolyte, 
where the positive and negative ions with concentrations c+ and c_ , respectively, have 
the same diffusivity D and valence charge number Z . 
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3.1. Bulk equations in the conventional ICEO model 

Neglecting chemical reactions in the bulk of the electrolyte, the ionic transport is 
governed by particle conservation through the continuity equation, 

where J± is the flux density of the two ionic species, respectively. Assuming a dilute 
electrolytic solution, the ion flux densities are governed by the Nernst-Planck equation, 

J±=-d(vc ± + ^c ± vA, (2) 



where the first term expresses ionic diffusion and the second term ionic electro- 
migration due to the electrostatic potential cf>. Here e is the elementary charge, T 
the absolute temperature and fee the Boltzmann constant. We note that due to the 
low fluid velocity v obtained in the ICEO systems under consideration, we can safely 
neglect the convective ion fluxes c±v throughout this paper, see Table [2 

The electrostatic potential <f> is determined by the charge density p c \ = Ze(c+—cJ) 
through Poisson's equation, 

V • (EfluidV^) = -/5cl, (3) 

where £fl u ;d is the fluid permittivity, which is assumed constant. The fluid velocity field 
v and pressure field p are governed the the continuity equation and the Navier-Stokes 
equation for incompressible fluids, 

V • v = 0, (4a) 
dv 



dt+ (v V)v 



- Vp + rjVv - p el Vcf), (46) 



where p m and rj are the fluid mass density and viscosity, respectively, both assumed 
constant. 



3.2. The artificial design field 7 used in the topology optimization model of ICEO 

To be able to apply the method of topology optimization, it is necessary to extend 
the conventional ICEO model with three additional terms, all functions of a position- 
dependent artificial design field 7(r). The design field varies continuously from zero 
to unity, where 7 = is the limit of a pure dielectric solid and 7 = 1 is the limit of a 
pure electrolytic fluid. The intermediate values of 7 represent a mixture of solid and 
fluid. 

The first additional term concerns the purely fluid dynamic part of our problem. 
Here, we follow Borrvall and Petersson [6] and model the dielectric solid as a porous 
medium giving rise to a Darcy friction force density —0(7)1?, where 0(7) may be 
regarded as a local inverse permeability, which we denote the Darcy friction. We let 
a (7) be a linear function of 7 of the form 01(7) = cx max (l —7), where a max = v/^poie i s 
the Darcy friction of the porous dielectric material assuming a characteristic pore size 
£ P orc- In the limit of a completely impenetrable solid the value of a max approaches 
infinity, which leads to a vanishing fluid velocity v. The modified Navier-Stokes 
equation extending to the entire domain, including the dielectric material, becomes 

dv 



ot + ^-^ 



= — Vp + rjV 2 v — p e lV(/> — a{^)v. (5) 



Gregersen et at : Topology and shape optimization of ICEO micropumps 



5 



The second additional term is specific to our problem. Since the Navier-Stokes 
equation is now extended to include also the porous dielectric medium, our model must 
prevent the unphysical penetration of the electrolytic ions into the solid. Following 
the approach of Kilic et al. |13j . where current densities are expressed as gradients 
of chemical potentials, J ex — V/i, we model the ion expulsion by adding an extra 
free energy term ^(7) to the chemical potential fi± — ±Zecf> + ksT ln(c± / cq) + ^(7) 
of the ions, where Co is the bulk ionic concentration for both ionic species. As above 
we let ^(7) be a linear function of 7 of the form ^(7) = K max (l — 7), where K max is 
the extra energy cost for an ion to enter a point containing a pure dielectric solid as 
compared to a pure electrolytic fluid. The value of K max is set to an appropriately 
high value to expel the ions efficiently from the porous material while still ensuring 
a smooth transition from dielectric solid to electrolytic fluid. The modified ion flux 
density becomes 

J± =-#(Vc ± + ^ c± V0+^c ± V K ( 7 )). (6) 

The third and final additional term is also specific to our problem. 
Electrostatically, the transition from the dielectric solid to the electrolytic fluid is 
described through the Poisson equation by a 7-dependent permittivity £(7). This 
modified permittivity varies continuously between the value Edici of the dielectric 
solid and £fl u id of the electrolytic fluid. As above, we would like to choose £(7) to 
be a linear function of 7. However, during our analysis using a metallic dielectric 
with £dici = 10 6 £o in an aqueous electrolyte with £fj u id = 78 £0 we found unphysical 
polarization phenomena in the electrolyte due to numerical rounding-off errors for 7 
near, but not equal to, unity. To overcome this problem wc ensured a more rapid 
convergence towards the value £fl u id by introducing a cut-off value 7 cu t — 0.98, a 
transition width A7 ~ 0.002, and the following expression for £(7), 

£(7) = £dici + (gfluid - £dici)|l - — 2 7 ^ tanh( 7c ^ 7 ) + 1 j- ( 7 ) 

F° r 7 ^ 7cut we obtain the linear relation £(7) = £di i + (£fluid — £dioi)7 7 while for 
7 ^ 7cut we have £(7) = £fl U id, see Fig. [UJb)-(c). For 7 sufficiently close to unity 
(and not only when 7 equals unity with numerical precision), this cut-off procedure 
ensures that the calculated topological break up of the dielectric solid indeed leads to 
several correctly polarized solids separated by pure electrolyte. The modified Poisson 
equation becomes 

V • [e(7)W] = -Pel- (8) 

Finally, we introduce the 7-dependent quantity, the so-called objective function 
$[7], to be optimized by the topology optimization routine: the flow rate in the x- 
direction perpendicular to the applied potential gradient. Due to incompressibility, 
the flow rate Q(x) is the same through cross-sections parallel to the yz-plane at any 
position x. Hence we can use the numerically more stable integrated flow rate as the 
objective function, 

$[7(r)] = / Q(x) dx — / v ■ n x dx dz, (9) 
Jo Jn 

where f2 is the entire geometric domain (including the design domain), and n x the 
unit vector in the x direction. 
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3. 3. Dimensionless form 

To prepare the numerical implementation, the governing equations are rewritten in 
dimensionless form, using the characteristic parameters of the system. In conventional 
ICEO systems the size a of the dielectric solid is the natural choice for the characteristic 
length scale £o, since the generated slip velocity at the solid surface is proportional 
to a. However, when employing topology optimization we have no prior knowledge 
of this length scale, and thus we choose it to be the fixed geometric half-length 
Co = H between the capacitor plates. Further characteristic parameters are the ionic 
concentration cq of the bulk electrolyte, and the thermal voltage (f>o — k&T / (Ze). 
The characteristic velocity uq is chosen as the Helmholtz-Smoluchowski slip velocity 
induced by the local electric field E = 4>o/Iq, an d finally the pressure scale is set by 
the characteristic microfluidic pressure scale po = r/uo/io. 

Although strictly applicable only to parallel-plate capacitors, the characteristic 
time To of the general system in chosen as the RC time of the double layer in terms 
of the Debye length Ad of the electrolyte [14] , 




T0= - Ad= «*'2(^- (10) 

Moreover, three characteristic numbers are connected to the 7-dependent terms 
in the governing equations: The characteristic free energy kq, the characteristic 
permittivity chosen as the bulk permittivity £fl u id, and the characteristic Darcy friction 
coefficient a a . In summary, 

l =H, 0o = -— , u = — , po = (11a) 

Ze r\ £ £ 

r = — — , uj = — , k = fc B T, ao = -jx- (Ho) 

L) TO tg 

The new dimensionless variables (denoted by a tilde) thus become 

(12a) 
(12b) 

In the following all variables and parameters are made dimensionless using these 
characteristic numbers and for convenience the tilde is henceforth omitted. 
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3.4- Linearized and reformulated equations 

To reduce the otherwise very time and memory consuming numerical simulations, we 
choose to linearize the equations. There are several nonlinearities to consider. 

By virtue of a low Reynolds number Re w 10 -6 , see Table El the nonlinear Navier- 
Stokes equation is replaced by the linear Stokes equation. Likewise, as mentioned in 
Sec. 13. H the low Peclet number Pe « 10~ 3 allows us to neglect the nonlinear ionic 
convection flux density c±v. This approximation implies the additional simplification 
that the electrodynamic problem is independent of the hydrodynamics. 

Finally, since ZeC/k^T < 0.5 the linear Debye-Hiickel approximation is valid, 
and we can utilize that the ionic concentrations only deviate slightly from the bulk 
equilibrium ionic concentration. The governing equations are reformulated in terms 
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of the average ion concentration c = (c + + c_)/2 and half the charge density 
p = (c_|_ — c_)/2. Thus, by expanding the fields to first order as c = 1 + 5c and 
p = + 5p, the resulting differential equation for p is decoupled from that of c. 
Introducing complex field notation, the applied external bias voltage is A<j)(t) = 
2VoCOs(u)t) = Rc[2Vo cxp(iwi)], yielding a corresponding response for the potential 
<j> and charge density p, with the complex amplitudes <f>(r) = 3>n(r) + i$/(r) and 
P(r) — Pfl(r) + iP/(r), respectively. The resulting governing equations for the 
electrodynamic problem is then 

V-[e(7)V$ fl ] =-^Pr, (13a) 

V-[e( 7 )V$ 7 ] = -^P/, (136) 

V • [V$ fi + VP R + P fl VK( 7 )] = -^P 7 , (13c) 

V • [V*j + VP/ + P/Vk( 7 )] = + ^Pr, (13d) 

where we have introduced the dimensionlcss thickness of the linear Dcbyc layer 
e = Ad/%- Given the electric potential $ and the charge density P, we solve for 
the time-averaged hydrodynamic fields (v) and (p), 

V • (u) = 0, (14a) 

= -V(p) + V 2 (t;) + (/ el ) - a(7)<«>, (146) 
where the time-averaged electric body force density (/ e i) is given by 

</el) = -^[i , JlV$fl + PjV$j]. (14c) 
5.5. Boundary conditions 

For symmetric dielectric solids we exploit the symmetry around z = and consider 
only the upper half (0 < z < 1) of the domain. As boundary condition on the driving 
electrode we set the amplitude Vq of the applied potential. Neglecting any electrode 
reactions taking place at the surface there is no net ion flux in the normal direction 
to the boundary with unit vector h. For the fluid velocity we set a no-slip condition, 
and thus at z — 1 we have 

- V , $/ - 0, (15a) 

n'[V$ fi + V^ + P fl V K ( 7 )]=0, (156) 

n- [V$/ + VP/ + P/Vk( 7 )] =0, (15c) 

(v) = 0. (15a") 

On the symmetry axis (z = 0) the potential and the charge density must be zero due 
to the anti-symmetry of the applied potential. Furthermore, there is no fluid flux in 
the normal direction and the shear stresses vanish. So at x — we have 

= = 0, Pr = Pi = 0, (16a) 

n-(v)=0, i-(tr)-n = 0, (166) 

where the dimensionless stress tensor is (<Tjfc) = —(p)5ik + (di(vk) + dk(vi)), and 
n and i are the normal and tangential unit vectors, respectively, where the latter 
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in 2D, contrary to 3D, is uniquely defined. On the remaining vertical boundaries 
(x — ±L/2£ ) periodic boundary conditions are applied to all the fields. 

Corresponding boundary conditions apply to the conventional ICEO model 
Eqs. (fTj)- (|4&|) , without the artificial design field but with a hard-wall dielectric solid. 
For the boundary between a dielectric solid and and electrolytic fluid the standard 
electrostatic conditions apply, moreover, there is no ion flux normal to the surface, 
and a no-slip condition is applied to the fluid velocity. 

4. Implementation and validation of numerical code 

4-1- Implementation and parameter values 

For our numerical analysis we use the commercial numerical finite-element modeling 
tool Comsol [IS] controlled by scripts written in Matlab [T5]. The mathematical 
method of topology optimization in microfluidic systems is based on the seminal paper 
by Borrvall and Petersson [6j, while the implementation containing the method of 
moving asymptotes by Svanberg [TTl HE] is taken from Olesen, Okkels and Bruus [7] . 

Due to the challenges discussed in Sec. 14.21 of resolving all length scales in the 
electrokinetic model, we have chosen to study small systems, 2H = 500 nm, with a 
relatively large Debye length, Ad = 20 nm. Our main goal is to provide a proof of 
concept for the use of topology optimization in electro-hydrodynamic systems, so the 
difficulties in fabricating actual devices on this sub-micrometric scale is not a major 
concern for us in the present work. A list of the parameter values chosen for our 
simulations is given in Table [T] 

For a typical topology optimization, as the one shown in Fig. (3Ja), approximately 
5400 FEM elements are involved. In each iteration loop of the topology optimization 



Table 1. Parameters used in the simulations of the topology optimization ICEO 
model and the conventional ICEO model. 



Parameter 


Symbol 


Dimensionless 


Physical 






value 


value 


Characteristic length 


to 


1.0 


250 nm 


Channel half-height 


H 


1.0 


250 nm 


Channel length 


L 


2.0 


500 nm 


Design domain half-height 


h 


0.8 


200 nm 


Design domain length 


I 


0.6 


150 nm 


Linear Debye length 


A D 


0.08 


20 nm 


Characteristic velocity 


Uq 


1.0 


1.7 x 10~ 3 m/s 


Characteristic potential 


<t>0 


1.0 


25 mV 


External potential amplitude 


v 


1.0 


25 mV 


External potential frequency 


u> 


6.28 


4 x 10 5 rad/s 


Bulk fluid permittivity 


Efluid 


1.0 


78 e 


Dielectric permittivity 


Ediel 


1.3 x 10 4 


10 6 e 


Bulk ionic concentration 


CO 


1.0 


0.23 mM 


Fluid viscosity 


V 


1.0 


10" 3 Pas 


Ionic diffusion constant 


D 


1.0 


2 x 10" 9 m 2 /s 


Ionic free energy in solid 


K 


3.0 


75 mV 


Maximum Darcy friction 




10 5 


2 x 10 16 Pas/m 2 
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routine three problems are solved: the electric problem , the hydrodynamic problem, 
and the adjunct problem for the sensitivity analysis, involving 4 x 10 4 , 2 x 10 4 , and 
7 x 10 4 degrees of freedom, respectively. On an Intel Core 2 Duo 2 GHz processer 
with 2 GB RAM the typical CPU time is several hours. 

4-2. Analytical and numerical validation by the conventional ICEO model 

We have validated our simulations in two steps. First, the conventional ICEO model 
not involving the design field j(r) is validated against the analytical result for the slip 
velocity at a simple dielectric cylinder in an infinite AC electric field given by Squires 
and Bazant [5] . Second, the design field model is compared to the conventional model. 
This two-step validation procedure is necessary because of the limited computer 
capacity. The involved length scales in the problem make a large number of mesh 
elements necessary for the numerical solution by the finite element method. Four 
different length scales appear in the gamma-dependent model for the problem of a 
cylinder placed mid between the parallel capacitor plates: The distance H from the 
center of the dielectric cylinder to the biased plates, the radius a of the cylinder, the 
Debye length Ad , and the length d over which the design field 7 changes from zero to 
unity. This last and smallest length-scale d in the problem is controlled directly be 
the numerical mesh size set up in the finite element method. It has to be significantly 
smaller than Ad to model the double layer dynamics correctly, so here a maximum 
value of the numerical mesh size is defined. 

The analytical solution of Squires and Bazant [3] is only strictly valid in the 
case of an infinitely thin Debye layer in an infinite electric field. So, to compare this 
model with the bounded numerical model the plate distance must be increased to 
minimize the influence on the effective slip velocity. Furthermore, it has been shown 
in a numerical study of Gregersen et al. [TO] that the Debye length Ad should be about 
a factor of 10 3 smaller than the cylinder radius a to approximate the solution for the 
infinitely thin Debye layer model. Including the demand of d being significantly smaller 
than Ad we end up with a length scale difference of at least I0 5 , which is difficult to 
resolve by the finite element mesh, even when mesh adaption is used. Consequently, 
we have checked that the slip velocity for the conventional model converges towards 
the analytical value when the ratio Ad /a decreases. Afterwards, we have compared 
the solutions for the conventional and gamma-dependent models in a smaller system 
with a ratio of Ad/o ~ 10 and found good agreement. 

4-3. Validation of the self- consistency of the topology optimization 

As an example of our validation of the self-consistency of the topology optimization, 
we study the dependence of the objective function Q — <&[u>, 7(w, r)] on the external 
driving frequency uj. As shown in Fig. [2ja)-(c) we have calculated the topology 
optimized dielectric structures jj = 7(0;., , r), j = a,b,c, for three increasing 
frequencies uj — uj a — 1.25, uj = ujb = 12.5, and uj — lj c = 62.5. In the following we let 
Qj{uj) denote the flow rate calculated at the frequency uj for a structure optimized at 
the frequency ujj. 

First, we note that Qj(ujj) decreases as the frequency increases above the 
characteristic frequency ujq — 1; Q a (uJ a ) — 2.95 x 10~ 3 , Qb(uJb) = 1-82 x ICP 3 , and 
Qc(uJ c ) — 0.55 x I0~ 3 . This phenomenon is a general aspect of ICEO systems, where 
the largest effect is expected to happen at uj — 2tt/tq = 6.28. 
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Figure 2. Validation of the self-consistency of the topology optimization for 
different driving frequencies w. (a) The streamline pattern (thick lines) for 
uj = ui a = 1-25 calculated using the design-field ICEO model with a porous 
dielectric medium (black and gray), the structure -y a of which has been found by 
topology optimization within the indicated rectangular design domain (straight 
lines). The flow rate for this converged solution structure is Q = 2.95 X 10~ 3 . (b) 
As panel (a) but with u) = = 12.5 and Q = 1.82 X 10 -3 . (c) As panel (a) but 
with u = uj c = 62.5 and Q = 0.55 X 10" 3 . (d) Flow rate Q versus frequency u 
for each of the three structures in panel (a), (b), and (c). Note that structure 7 a 
indeed yields the highest flow rate Q for w = w a , structure ^ maximizes Q for 
u) = ii)\y. and structure -y c maximizes Q for u) = u) c . 

Second, and most significant, we see in Fig. [2fd) that structure 7 a is indeed the 
optimal structure for u = uj a since Q a (w a ) > Qb(w a ), Q c {u a )- Likewise, 7;, is optimal 
for u> — LUb, and j c is optimal for to = lu c . 

We have gained confidence in the self-consistency of our topology optimization 
routine by carrying out a number of tests like the one in the example above. 

5. Results 

5.1. Topology optimization 

For each choice of parameters the topology optimization routine converges to a specific 
distribution of dielectric solid given by 7(1"). As a starting point for the investigation of 
the optimization results we used the parameters listed in Table [T] As discussed above, 
the geometric dimensions are chosen as large as possible within the computational 
limitations: the Debye length is set to Ad = 20 nm and the distance between the 
capacitor plates to 2H = 500 nm. The external bias voltage is of the order of 
the thermal voltage Vq = 25 mV to ensure the validity of the linear Debye-Huckel 
approximation. We let the bulk fluid consist of water containing small ions, such as 
dissolved KC1, with a concentration c = 0.23 mM set by the chosen Debye length. 
The dielectric material permittivity is set to e^ici = 10 6 £0 m order to mimic the 
characteristics of a metal. The artificial parameters k and a max are chosen on a pure 
computational basis, where they have to mimic the real physics in the limits of fluid 
and solid, but also support the optimization routine when the phases are mixed. 

Throughout our simulations we have primarily varied the applied frequency uj 
and the size I X 2/i of the design domain. In Fig. [2] we have shown examples of large 
design domains with Ix h = 2.0 x 0.8 covering 80% of the entire domain and frequency 
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Figure 3. (a) The streamline pattern (thick lines) calculated for u> = 6.28 
using the design-field ICEO model with a porous dielectric medium (black and 
gray), the structure of which has been find by topology optimization within the 
indicated rectangular design domain (thin lines). The flow rate for this converged 
solution structure is Q = 2.99 X 10~ 3 . (b) The streamline pattern (full lines) 
calculated using the conventional ICEO model with a hard- walled dielectric solid 
(black). The shape of the dielectric solid is the 0.95-contour of the 7- field taken 
from the topology-optimized structure shown in panel (a). The flow rate is 
Q = Q* = 1.88 X 10" s . (c) and (d) Color plot of the charge density p(r) 
corresponding to panel (a) and (b), respectively. See Table[T]for parameter values. 



sweeps over three orders of magnitude. However, in the following we fix the frequency 
to be w = 27t/to = 6.28, where the ICEO response is maximum. Moreover we focus 
on a smaller design domain I x h — 0.6 x 0.8 to obtain better spatial resolution for 
the given amount of computer memory and to avoid getting stuck in local optima. It 
should be stressed that the size of the design domain has a large effect on the specific 
form and size of the dielectric islands produced by the topology optimization. Also, 
it is important if the design domain is allowed to connect to the capacitor plates or 
not, see the remarks in Sec. [5] 

The converged solution found by topology optimization under these conditions 
is shown in Fig. 02a). The shape of the porous dielectric material is shown together 



Table 2. The value of characteristic physical quantities calculated in the topology 
optimization ICEO model corresponding to Fig. \3\ 



Quantity 


Symbol 


Dimcnsionless 


Physical 






value 


value 


Gap between dielectric pieces 


c gap 


0.4 


100 nm 


Velocity in the gap 


Ugap 


0.016 u 


28 /xm/s 


Largest zeta potential 


Cm ax 


0.5^o 


12.5 mV 


Reynolds number Re 


Pm^gap^gap/^7 


2.8 x 10~ 6 




Pcclct number Pe 


^gap^gap/-^ 


1.4 x 10~ 3 




Debye-Hiickel number Hii 


eCmax/(4fc B T) 


0.13 
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with a streamline plot of equidistant contours of the flow rate. We notice that many 
stream lines extend all the way through the domain from left to right indicating that 
a horizontal flow parallel to the ir-axis is indeed established. The resulting flow rate 
is Q = 2.99 x 10~ 3 . The ICEO flow of this solution, based on the design-field model, 
is validated by transferring the geometrical shape of the porous dielectric medium 
into a conventional ICEO model with a hard-walled dielectric not depending on the 
design field. In the latter model the sharp interface between the dielectric solid and 
the electrolyte is defined by the 0.95-contour of the topology optimized design field 
j(r). The resulting geometry and streamline pattern of the conventional ICEO model 
is shown in Fig. 0Jb). The flow rate is now found to be Q = Q* = 1.88 x 1CT 3 . 
There is a close resemblance between the results of two models both qualitatively and 
quantitatively. It is noticed how the number and positions of the induced flow rolls 
match well, and also the absolute values of the induced horizontal flow rates differs 
only by 37%. 

Based on the simulation we can now justify the linearization of our model. The 
largest velocity u gap is found in the gap of width £ gap between the two satellite 
pieces and the central piece. As listed in Table [2] the resulting Reynolds number 
is Re = 2.8 x 10~ 7 , the Peclet number is Pe — 1.4 x 1CP 3 , while the Debye-Hiickel 
number is Hii = 0.13. 

5.2. Comparison to simple shapes 

We evaluate our result for the optimized flow rate by comparing it to those obtained 
for more basic, simply connected, dielectric shapes, such as triangles and perturbed 
circles previously studied in the literature as flow inducing objects both analytically 
and experimentally [31 01 [5] ■ For comparison, such simpler shapes have been fitted 
into the same design domain as used in the topology optimization Fig. [3fa), and the 
conventional ICEO model without the design field was solved for the same parameter 
set. In Fig. Eta) the resulting flow for a triangle with straight faces and rounded 
corners is shown. The height b of the face perpendicular to the symmetry line was 
varied within the height of the design domain < b < 0.8, and the height b = 0.32 
generating the largest flow in the ir-direction results in a flow rate of Q = 0.22 x 10 -3 , 
which is eight times smaller than the topology optimized result. In Fig. 0Jb) the 




Figure 4. (a) The streamline pattern (thick lines) for a simple triangular 
reference structure calculated for oj = 6.28 using the conventional ICEO model 
with a hard-walled dielectric solid (black). The height b = 0.32 of the triangle is 
chosen to give the largest flow rate for a fixed base line given by the rectangular 
design domain of Fig. |3{a). The flow rate is Q = 0.22 X 10~ 3 . (b) The same 
as panel (a) except the geometry of the dielectric solid is given by the perturbed 
circle r(0) = 0.24[1 + 0.5 cos(30)] . The flow rate is Q = 0.46 X 10" 3 . 
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induced flow around a perturbed cylinder with radius r(0) = 0.24 [l + 0.5cos(3#)] is 
depicted. Again the shape has been fitted within the allowed design domain. The 
resulting flow rate Q = 0.46 x 10~ 3 is higher than for the triangle but still a factor 
of four slower than the optimized result. It is clearly advantageous to change the 
topology of the dielectric solid from simply to multiply connected. 

For the topology optimized shape in Fig. [3][a) it is noticed that only a small 
amount of flow is passing between the two closely placed dielectric islands in the upper 
left corner of the design domain. To investigate the importance of this separation, the 
gap between the islands was filled out with dielectric material and the flow calculated. 
It turns out that this topological change only lowered the flow rate slightly (15%) to 
a value of Q = 1.59 x 10~ 3 . Thus, the important topology of the dielectric solid in 
the top-half domain is the appearance of one center island crossing the antisymmetry 
line and one satellite island near the tip of the center island. 

5.3. Shape optimization 

The topology optimized solutions are found based on the extended ICEO model 
involving the artificial design field j(r). To avoid the artificial design field it is 
desirable to validate and further investigate the obtained topology optimized results 
by the physically more correct conventional ICEO model involving hard-walled solid 
dielectrics. We therefore extend the reasoning behind the comparison of the two 
models shown in Fig. [3] and apply a more stringent shape optimization to the various 
topologies presented above. With this approach we are gaining further understanding 
of the specific shapes comprising the overall topology of the dielectric solid. Moremore, 
it is possible to point out simpler shapes, which are easier to fabricate, but still perform 
well. 

In shape optimization the goal is to optimize the objective function $, which 
depends on the position and shape of the boundary between the dielectric solid and 
the electrolytic fluid. This boundary is given by a line interpolation through a small 
number of points on the boundary. These control points are in turn given by N design 
variables g — (51,32, ■ ■ ■ > <?iv); so the objective function of Eq. ([9]) depending on the 
design field "f(r) is now written as <&[«/] depending on the design variables g, 

= / v ■ n x dx dz. (17) 
Jn 

To carry out the shape optimization we use a direct bounded Nelder-Mead simplex 
method [20] implemented in Matlab [2lJ [22] . This robust method finds the optimal 
point <jr op t in the .ZV-dimensional design variable space by initially creating a simplex 
in this space, e.g. a iV-dimensional polyhedron spanned by N + 1 points, one of which 
is the initial guess. The simplex then iteratively moves towards the optimal point 
by updating one of the N + 1 points at the time. During the iteration, the simplex 
can expand along favorable directions, shrink towards the best point, or have its 
worst point replaced with the point obtained by reflecting the worst point through 
the centroid of the remaining N points. The iteration terminates once the extension 
of the simplex is below a given tolerance. We note that unlike topology optimization, 
the simplex method relies only on values of the objective function and not on the 
sensitivity d^/dg [23] . 

First, we perform shape optimization on a right-angled triangle corresponding to 
the one shown in Fig. SJa). Due to the translation invariance in the cc-direction, we 
fix the first basepoint of the triangle (xi,0) to the right end of the simulation domain, 
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Figure 5. (a) The streamline pattern (thick lines) for the shape-optimized right- 
angled triangle fixed at the symmetry line 2 = calculated for w = 6.28 using 
the conventional ICEO model with a hard- walled dielectric solid (black). In 
the full domain this is a triangle symmetric around 2 = 0. The flow rate is 
Q = 0.32 X 10 — 3 . (b) As in panel (a) but without constraining the triangle to be 
right-angled. In the full domain the shape is foursided polygon symmetric around 
2 = 0. The flow rate is Q = 0.76 X 10 — 3 . Note that all sharp corners of the 
polygons have been rounded by circular arcs of radius 0.01. 



while the second point (2:2,0) can move freely along the baseline, in contrast to the 
original rectangular design. To ensure a right-angled triangle only the z-coordinatc 
of the top point (£2, £2) may move freely. In this case the design variable becomes 
the two-component vector g = (2:2, #2)- The optimal right-angled triangle is shown 
in Fig. OUa). The flow rate is Q = 0.32 x 10~ 3 or 1.5 times larger than that of the 
original right-angled triangle confined to the design domain. 

If we do not constrain the triangle to be right-angled, we instead optimize a 
polygon shape spanned by three corner points in the upper half of the electrolytic 
capacitor. So, due to the symmetry of the problem, we are in fact searching for the 
most optimal, symmetric foursided polygon. The three corner points are now given 
as (xj.,0), (22,0), and (23,2:3), and again due to translation invariance, it results in 
a three-component design variable g — (22,23,23). The resulting shape-optimized 
polygon is shown in Fig. ODJb). The flow rate is Q — 0.76 x 10~ 3 , which is 3.5 times 
larger than that of the original right-angled triangle confined to the design domain 
and 2.4 times better than that of the best right-angled triangle. However, this flow 
rate is still a factor of 0.4 lower than the topology optimized results. 

To be able to shape optimize the more complex shapes of Fig. [3] we have employed 
two methods to obtain a suitable set of design variables. The first method, the radial 
method, is illustrated in Fig. [6] The boundary of a given dielectric solid is defined 
through a cubic interpolation line through N control points (xi,Zi),i = 1,2, . . . ,N, 
which are parameterized in terms of two co-ordinates (x c , z c ) of a center point, two 
global scale factors A and B, N lengths rj, and N fixed angles Qi distributed in the 
interval from to 2ir, 



In this case the design variable becomes g = (x c , z c , r\,r%, . . . , r/v, A, B). 

The second parametrization method involves a decomposition into harmonic 
components. As before we define a central point (x c , z c ) surrounded by N control 
points. However now, the distances are decomposed into M harmonic components 
given by 



(Xi,Zi) = (2 c ,z c ) + r, i (^4cos6>i, Bsin^). 



(18) 




n=l 



(19) 
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Figure 6. Illustration of the parametrization, Eq. JTH}, of the boundary 
of a dielectric solid with a complex shape. The polar representation is 
shown for point i = 7. The shape consist of five harmonic components 
represented by Eq. H19H with the design- variables x c = —0.1312, z c = 
0.7176, r = 0.1403, A> = {0.2501,0.0151,0.0062,0.2103,0.2313}, ip t = 
{-1.7508, -2.2526, 0.4173, 0.1172, -0.2419}. 

where r$ is an overall scale parameter and ip n is a phase shift. In this case the design 
variable becomes g = (x c , z c , r , Ai,A 2 , A M , ¥>i, ¥>2, ■ ■ ■ , <Pm)- 

5.4- Comparing topology optimization and shape optimization 

When shape-optimizing a geometry similar to the one found by topology optimization, 
we let the geometry consist of two pieces: (i) an elliptic island centered on the 
symmetry-axis and fixed to the right side of design domain, and (ii) an island with a 
complex shape to be placed anywhere inside the design domain, but not overlapping 
with the elliptic island. For the ellipse we only need to specify the major axis A and 
the minor axis B, so these two design parameters add to the design variable listed 
above for either the radial model or the harmonic decomposition model. To be able 
to compare with the topology optimized solution the dielectric solid is restricted to 
the design domain. 

The result of this two-piece shape optimization is shown in Fig. [7] Compared to 
the simply connected topologies, the two-piece shape-optimized systems yields much 
improved flow rates. For the shape optimization involving the radial method with 16 
directional angles and A = B for the complex piece, the flow rate is Q — 1.92 x 10 -3 , 
Fig. Ha), which is 2.5 times larger than that of the shape-optimized foursided 
symmetric polygon. The harmonic decomposition method, Fig. E^b), yields a flow 
rate of Q = 1.52 x 1CP 3 or 2.0 times larger than that of the polygon. 

All the results for the obtained flow rates are summarized in Table [3l It is seen 
that two-piece shape optimized systems performs as good as the topology optimized 
system, when analyzed using the conventional ICEO model without the artificial 
design field. We also note by comparing Figs.[3]and[7]that the resulting geometry found 
using either topology optimization or shape optimization is essentially the same. The 
central island of the dielectric solid is a thin structure perpendicular to the symmetry 
axis and covering approximately 60% of the channel width. The satellite island of 
complex shape is situated near the tip of the central island. It has two peaks pointing 
towards the central island that seem to suspend a flow roll which guides the ICEO 
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Figure 7. Shape-optimized dielectrics with a topology corresponding to the 
topology-optimized shapes of Fig. [3] (a) The streamline pattern (thick lines) for 
a two-piece geometry calculated using the conventional ICEO model. The shape 
of the hard-walled dielectric solid (black) is found by shape optimization using 
the radial method Eq. H18H with N = 16 directional angles. The flow rate is 
Q = 1.92 X 10~ 3 . (b) The same as panel (a) except the geometry of the dielectric 
solid is by shape optimization using the harmonic decomposition method Eq. H19I I 
with M = 5 modes. The flow rate is Q = 1.52 X 10~ 3 . (c) and (d) Color plot of 
the charge density p(r) corresponding to panel (a) and (b), respectively. 



flow through the gap between the two islands. 



6. Concluding remarks 

The main result of this work is the establishment of the topology optimization method 
for ICEO models extended with the design field j(r). In contrast to the conventional 
ICEO model with its sharply defined, impenetrable dielectric solids, the design field 
ensures a continuous transition between the porous dielectric solid and the electrolytic 



Table 3. Overview of the resulting flow rates Q relative to the topology optimized 
value Q* = 1.88 X 10 3 , see Fig. (3{b), for the various geometries analyzed in the 
conventional ICEO model. The methods by which the geometries have been 
determined are listed. 



Shape 


Method 


Flow rate 






Q/Q* 


Triangle with optimal height, Fig. Ufa) 


Shape optimization 


0.12 


Perturbed cylinder, Fig. 0|b) 


Fixed shape 


0.24 


Optimized triangle, Fig. OJa) 


Shape optimization 


0.17 


Optimized foursided polygon, Fig. [5jb) 


Shape optimization 


0.40 


Topology optimized result, Fig. [3£b) 


Topology optimization 


1.00 


Harmonic decomposition and ellipse, Fig. [7]^ a) 


Shape optimization 


0.81 


Radial varying points and ellipse, Fig. [7jb) 


Shape optimization 


1.02 



Gregersen et ai: Topology and shape optimization of ICEO micropumps 17 

fluid, which allows for an efficient gradient-based optimization of the problem. By 
concrete examples we have shown how the use of topology optimization has led to non- 
trivial system geometries with a flow rate increase of nearly one order of magnitude. 

However, there exist many local optima, and we cannot be sure that the converged 
solution is the global optimum. The resulting shapes and the generated flow rates 
depend on the initial condition for the artificial 7-field. Generally, the initial condition 
used throughout this paper, 7 = 0.99 in the entire design domain, leads to the most 
optimal results compared to other initial conditions. This initial value corresponds 
to a very weak change from the electrolytic capacitor completely void of dielectric 
solid. In contrast, if we let 7 — 0.01 corresponding to almost pure dielectric 
material in the entire design region, the resulting shapes are less optimal, i.e. the 
topology optimization routine is more likely to be caught caught in a local optimum. 
Furthermore, the resulting shapes turns out to be mesh-dependent as well. So, 
we cannot conclude much about global optima. Instead, we can use the topology 
optimized shapes as inspiration to improve existing designs. For this purpose shape 
optimization turns out to be a powerful tool. We have shown in this work how shape 
optimization can be used efficiently to refine the shape of the individual pieces of the 
dielectric solid once its topology has been identified by topology optimization. 

For all three additional 7-dependent fields 01(7), ^(7), and £(7) we have used 
(nearly) linear functions. In many previous applications of topology optimization 
non-linear functions have successfully been used to find global optima by gradually 
changing the non-linearity into strict linearity during the iterative procedure [SI El El 
112) . However, we did not improve our search for a global optimum by employing such 
schemes, and simply applied the (nearly) linear functions during the entire iteration 
process. 

The limited size of the design domain is in some cases restricting the free formation 
of the optimized structures. This may be avoided by enlarging the design domain. 
However, starting a topology optimization in a very large domain gives a huge amount 
of degrees of freedoms, and the routine is easily caught in local minima. These local 
minima often yield results not as optimal as those obtained for the smaller design 
boxes. A solution could be to increase the design domain during the optimization 
iteration procedure. It should be noted that increasing the box all the way up to 
the capacitor plates results in solution shapes, where some of the dielectric material 
is attached to the electrode in order to extend the electrode into the capacitor and 
thereby maximize the electric field locally. This may be a desirable feature for some 
purposes. In this work we have deliberately avoided such solutions by keeping the 
edges of the design domain from the capacitor plates. 

Throughout the paper we have only presented results obtained for dielectric 
solids shapes forced to be symmetric around the center plane z = 0. However, we 
have performed both topology optimization and shape optimization of systems not 
restricted to this symmetry. In general we find that the symmetric shapes always 
are good candidates for the optimal design. It cannot be excluded, though, that in 
some cases a spontaneous symmetry breaking occurs similar to the asymmetric S-turn 
channel studied in Ref. [7]. 

By studying the optimized shapes of the dielectric solids, we have noted that 
pointed features often occurs, such as those clearly seen on the dielectric satellite 
island in Fig. EJb). The reason for these to appear seems to be that the pointed 
regions of the dielectric surfaces can support large gradients in the electric potential 
and associated with this also with large charge densities. As a result large electric 
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body forces act on the electrolyte in these regions. At the same time the surface 
between the pointed features curve inward which lowers the viscous flow resistance 
due to the no-slip boundary condition. This effect is similar to that obtained by 
creating electrode arrays of different heights in AC electro-osmosis [MJ [25] . 

Another noteworthy aspect of the topology optimized structures is that the 
appearance of dielectric satellite islands seem to break up flow rolls that would 
otherwise be present and not contribute to the flow rate. This leads to a larger 
net flow rate, as can be seen be comparing Figs. |4] and [7l 

Throughout the paper have treated the design field 7 as an artificial field. 
However, the design-field model could perhaps achieve physical applications to systems 
containing ion exchange membranes, as briefly mentioned in Sec.[T] Such membranes 
are indeed porous structures permeated by an electrolyte. 

In conclusion, our analysis points out the great potential for improving actual 
ICEO-based devices by changing simply connected topologies and simple shapes of 
the dielectric solids, into multiply connected topologies of complex shapes. 
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